#include "Pyramid.h"
#include "helper.h"
#include "string.h"
#include "include/stx/btree_multimap.h"


#include <cstdlib> // for exit function
#include <vector>


typedef double key_type;
typedef int value_type;

struct traits_nodebug
{
    static const bool       selfverify = false;
    static const bool       debug = false;

    static const int        leafslots = 8;
    static const int        innerslots = 8;
};

typedef stx::btree_multimap<key_type, value_type,std::less<key_type>, struct traits_nodebug> btree_type;

btree_type bt;

inline std::string trimStr(std::string& str)
{
  str.erase(0, str.find_first_not_of(" \t\r\n"));
  str.erase(str.find_last_not_of(" \t\r\n")+1);

  return str;
}
//A function to set number of dimensions
void Pyramid::setDimension(int d){
    D=d;
}

void Pyramid::setdata(double ** d,int tr, int dimension ){

    D=dimension;
    totalrecords=tr;
    data=d;
    pointID = new int[tr];
    pv = new pyrval[tr];



}

void Pyramid::savemedian()
{
	//cout << "saving median "<<endl;
	for(int i=0; i<D; i++)
	{

		outFile << am.getmedian(i) << endl;
		//cout << am.getmedian(i) << endl;
	}

}

void Pyramid::openOutFile(char * fileName)
{
    outFile.open(fileName);
}

void Pyramid::closeOutFile()
{
	outFile.close();
}

void Pyramid::openInFile(char * fileName)
{
    inFile.open(fileName, ios::in);
    inFile.precision(3);
}

void Pyramid::closeInFile()
{
	inFile.close();
}

void Pyramid::saveblankmedian()
{

	for(int i=0; i<D; i++)
	{
		outFile << "0.5" << endl;
	}

}

void Pyramid::loadmedian()
{
	am.initialize(totalrecords, D, data);


 //cout << "Loading median " << endl;
 string line;
 int lineI = -1;

 if (inFile.is_open())
 {

   while ( inFile.good() )
   {
     getline(inFile,line);
     line = trimStr(line);
     if(line.length() != 0)
     {
		 lineI = lineI + 1; //first line is 0

		 char* lineStr = new char[line.length()+1]; //+1 for auto added null terminator
		 strcpy(lineStr, line.c_str());

		 am.set(lineI, strtod(lineStr, NULL));
		 //cout << lineStr << endl;
	  }
   }
 }
 else
 {
   cout << "Unable to open file" << endl;
 }
}


void Pyramid:: buildSequentialScanData()
{
	if(data== NULL)
	{
		cout<<"Data not loaded \n";
		return;
	}

	for(int i=0; i< totalrecords; i++)
	{
		pointID[i] =  data[i][0]; //save off data point id

	}
}
void Pyramid:: buildindexAM()
{

    if(data== NULL)
    {
        cout<<"Data not loaded \n";
        return;
    }

    double start = startTimer();
    double *array_data = new double[totalrecords*D];
    double **cvdata = new double*[totalrecords];
    for (int i = 0; i < totalrecords; ++i)
    cvdata[i] = array_data + D*i;

    for(int i=0; i< totalrecords; i++){
    	pointID[i] =  data[i][0]; //save off data point id
        for(int j=0; j<D; j++)
        {
            cvdata[i][j] = data[i][j+1];
        }
    }

    algam(totalrecords, D, cvdata);

    for(int i=0; i< totalrecords; i++)
    {
        pv[i] = pyramidvalue_am(cvdata[i]);
    }
    double stop = stopTimer();
    indexTime = stop - start;
    outFile << "Time to build index: " << indexTime << " Seconds " << endl;
}


void Pyramid::buildindexM()
{
   if(data== NULL)
    {
        cout<<"Data not loaded \n";
        return;
    }
    double start = startTimer();
    double *array_data = new double[totalrecords*D];
    double **cvdata = new double*[totalrecords];
    for (int i = 0; i < totalrecords; ++i)
    cvdata[i] = array_data + D*i;

    for(int i=0; i< totalrecords; i++){
    	pointID[i] =  data[i][0]; //save off data point id
        for(int j=0; j<D; j++)
        {
            cvdata[i][j] = data[i][j+1];
        }
    }

    findmedian(totalrecords, D, cvdata);

    for(int i=0; i< totalrecords; i++)
    {
        pv[i] = pyramidvalue_m(cvdata[i]);
    }
    double stop = stopTimer();
    indexTime = stop - start;
    outFile << "Time to build index: " << indexTime << " Seconds " << endl;
}

void Pyramid:: buildindexHM()
{
    if(data== NULL)
    {
        cout<<"Data not loaded \n";
        return;
    }

    double start = startTimer();
    double *array_data = new double[totalrecords*D];
    double **cvdata = new double*[totalrecords];
    for (int i = 0; i < totalrecords; ++i)
    cvdata[i] = array_data + D*i;


    for(int i=0; i< totalrecords; i++){
    	pointID[i] =  data[i][0]; //save off data point id
        for(int j=0; j<D; j++)
        {
            cvdata[i][j] = data[i][j+1];
        }
    }
   
    hbam(totalrecords, D, cvdata);
   
    for(int i=0; i< totalrecords; i++)
    {
        pv[i] = pyramidvalue_am(cvdata[i]);
        //if(i==50)
        	//cout << "Pyramid value for " << pointID[i] << " is " << pv[i].pyramid + pv[i].height << endl;
    }
    double stop = stopTimer();
    indexTime = stop - start;
    outFile << "Time to build index: " << indexTime << " Seconds " << endl;

}
void Pyramid::buildindex(){
    if(data == NULL){
        cout<<"Data not loaded \n";
        return;
    }

    double start = startTimer();

    for(int i=0; i< totalrecords; i++){
        double tempdata[D];

        pointID[i] =  data[i][0]; //save off data point id

        for(int j=0; j< D; j++){
            tempdata[j]= data[i][j+1]; // skip over data point id
        }
        pv[i] = pyramidvalue(tempdata);

    }
    double stop = stopTimer();
    indexTime = stop - start;
    outFile << "Time to build index: " << indexTime << " Seconds " << endl;
}

void Pyramid::buildtree(){

	double start = startTimer();
	for(int i=0;i <totalrecords; i++)
	{

		bt.insert(pv[i].pyramid + pv[i].height, i );
	}

	double stop = stopTimer();
	treeTime = stop - start;
	outFile << "Time to build Tree: "  << treeTime << " Seconds" << endl;
	outFile << "Total build time: "  << (treeTime + indexTime) << " Seconds" << endl;

	outFile << "Size of tree: " <<bt.get_stats().levels << " levels, " << bt.get_stats().nodes() << " nodes."<< endl;
	outFile << "Tree segments:\n";
	bt. print_tree_segments(outFile);



}

void Pyramid::loadtree(){

	bt.restore(inFile);
	//printTree();
    for(int i=0; i< totalrecords; i++){
    	pointID[i] =  data[i][0]; //save off data point id
    }

}

void Pyramid::savetree(){
	//printTree();
	bt.dump(outFile);
}


void Pyramid::printTree()
{
	for(btree_type::iterator bi = bt.begin(); bi != bt.end(); bi++)
	{
		cout << "Key: " << bi.key() << " data: " <<bi.data() << endl;
	}
}
void Pyramid::searchquerypoints(double **d,int tr)
{

    pyrval pv;

    double * temp = new double[D];

    for(int i=0; i< tr; i++)
    {
    	double start = startTimer();
    	for(int j=0;j<D;j++)
    	{
    		temp[j] = d[i][j+1];
    	}
        pv=pyramidvalue(temp);

        outFile << "Point(s) found: ";
    	stx::CountQueryNodes("Start");
    	btree_type::iterator bi = bt.lower_bound(pv.pyramid + pv.height);
        while(bi != bt.end() && bi.key() <= pv.pyramid + pv.height){

        	outFile << pointID[bi.data()] << " ";

        	bi++;
        }
        outFile << endl;

        outFile<< "Total nodes searched : " <<stx::CountQueryNodes("Get")<<endl;
    	double stop = stopTimer();
    	outFile <<"Query Time: " << (stop - start)<< " Seconds\n";
    	outFile << endl;

    }



}

void Pyramid::searchquerypointsUAppM(double **d,int tr)
{

    pyrval pv;

    double * temp = new double[D];

    for(int i=0; i< tr; i++)
    {
    	double start = startTimer();
    	for(int j=0;j<D;j++)
    	{
    		temp[j] = d[i][j+1];
    	}
        pv=pyramidvalue_am(temp);

        outFile << "point(s) found: ";
    	stx::CountQueryNodes("Start");
    	btree_type::iterator bi = bt.lower_bound(pv.pyramid + pv.height);
        while(bi != bt.end() && bi.key() <= pv.pyramid + pv.height){

        	outFile << pointID[bi.data()] << " ";

        	bi++;
        }
        outFile << endl;

        outFile<< "Total nodes searched : " <<stx::CountQueryNodes("Get")<<endl;
        double stop = stopTimer();
        outFile <<"Query Time: " << (stop - start)<< " Seconds\n";
        outFile << endl;
    }

}

void Pyramid::searchqueryrangesUAppM(double* min, double* max, int tr)
{
	stx::CountQueryNodes("Start");
    double qr[D][2];
    double lower;
    double upper;

    for(int line=0;line<tr;line++)
	{
    	double start = startTimer();
    	stx::CountQueryNodes("Start");
    	outFile << "Range Query: " << endl;
		for(int i=0; i< D; i++)
		{
			qr[i][0] = min[(line*D) + i];
			outFile << min[(line*D) + i] << ",";
			qr[i][1] = max[(line*D) + i];
			outFile << max[(line*D) + i] << "\n";

		}


	    //Use the query range to perform KNN search
	    Pyramid::queryrange *q = getqueryrange_am(qr);
	    //printTree();
	    outFile << "Points in range: ";
	    for(int i=0; i<2*D;i++)
	    {

	    	if(q[i].pyramid != -1)
	    	{
	    		lower = q[i].pyramid + q[i].h_low;
	    		upper = q[i].pyramid + q[i].h_high;
	    		//cout << "pyramid " << i << " hlow " << q[i].h_low<< " h_high " << q[i].h_high << endl;
	    		btree_type::iterator bi = bt.lower_bound(lower);
	    		//cout << "bi.data(): " << bi.data() << " bt.end() " << endl;
	    	    while(bi != bt.end() && bi.key() <= upper)
	    	    {
	    	    	//cout << "bi.data(): " << bi.data() << endl;
	    	    	if(inside(qr, bi.data()))
	    	    	{
	    	    		outFile << pointID[bi.data()] << " ";
	    	    		//cout << "point in range: " << pointID[bi.data()] << endl;
	    	    	}
	    	    	bi++;
	    	    }
	    	}
	    }

	    outFile << "\n";

	    outFile<< "Total nodes searched : " <<stx::CountQueryNodes("Get")<<endl;
	    double stop = stopTimer();
	    outFile <<"Query Time: " << (stop - start)<< " Seconds\n";
	    outFile << endl;
	}


}

void Pyramid::searchqueryKnnUAppM(double **d, int tr,int k)
{
	for(int line =0;line<tr;line++)
	{
		double start = startTimer();
			 stx::CountQueryNodes("Start");
			 double dataPoint[D];
			 for(int i=0;i<D;i++)
			 {
				 dataPoint[i] = d[line][i+1]; //strip off point id
			 }

			 Pyramid::pyrval pv = pyramidvalue(dataPoint);
			 double qr[D][2];
			 double qw[D][2];
			 double ** convertedQr;
			 double tempLargestValue = 0;
			 int 	tempLargestIndex;
			 vector<int> A;
			 double dMax = 0;

			 searchLeft(pv, A, dataPoint, k, &dMax);
			 searchRight(pv, A, dataPoint, k, &dMax);

			 //did we find k points in first pyrimid?
			 if( A.size() < k)
			 {
				 //make qr entire data space
				for(int i=0;i<D;i++)
				{
					qr[i][0] = 0.0;
					qr[i][1] = 1.0;
				}
			 }

			 for(int j=0; j<2*D; j++)
			 {
				 if(j != pv.pyramid)
				 {
					 //shrink qr if we have found k candidates
					 if( A.size() >= k)
					 {
						 for(int k=0;k<D;k++)
						 {

							 //calculate min in the kth dimension
							 if(dataPoint[k] - dMax < 0)
							 {
								 qr[k][0] = 0;
							 }
							 else
							 {
								 qr[k][0] = dataPoint[k] - dMax;
							 }

							 //calculate max in the kth dimension
							 if(dataPoint[k] + dMax > 1)
							 {
								 qr[k][1] = 1;
							 }
							 else
							 {
								 qr[k][1] = dataPoint[k] + dMax;
							 }
						 }
					 }


					 //find intersection of qr with j
					 convertedQr = am.convertqueryrange(qr);
					 double hlow, hhigh;
					 for(int i=0; i<D; i++)
					 {
				 		 qw[i][0] = convertedQr[i][0] - 0.5;
						 qw[i][1] = convertedQr[i][1] - 0.5;
					 }
					 if(intersect(qw,j,&hlow,&hhigh))
					 {
						 //update A and dMax
						double lower = j + hlow;
						double upper = j + hhigh;
						btree_type::iterator bi = bt.lower_bound(lower);
						while(bi != bt.end() && bi.key() <= upper)
						{
							//check if bi.data() needs to be added to A
							if(A.size() < k)
							{
								 //add point to A then check if we need to update dMax
								 A.push_back(bi.data());

								 if(euclideandistance(&data[bi.data()][1], dataPoint) > dMax )
								 {
									 dMax = euclideandistance(&data[bi.data()][1], dataPoint);
								 }
							}
							else
							{
								//if distance between bi  and q < dMax update A and dMax
								 if(euclideandistance(&data[bi.data()][1], dataPoint) < dMax && !in(A, bi.data()))
								 {
									 //update A
									 tempLargestValue = 0;
									 for(int i=0;i<A.size(); i++)
									 {
										 if(euclideandistance(&data[A[i]][1], dataPoint) > tempLargestValue)
										 {
											 tempLargestValue = euclideandistance(&data[A[i]][1], dataPoint);
											 tempLargestIndex = i;
										 }
									 }
									 A[tempLargestIndex] = bi.data();

									 //update dMax
									 tempLargestValue = 0;
									 for(int i=0;i<A.size(); i++)
									 {
										 if(euclideandistance(&data[A[i]][1], dataPoint) > tempLargestValue)
										 {
											 tempLargestValue = euclideandistance(&data[A[i]][1], dataPoint);
										 }
									 }
									 dMax = tempLargestValue;
								 }
							}

							bi++;
						}
					 }
				 }
			 }


			 outFile << "k-NN Query " << endl;
			 double stop = stopTimer();

			 //print out query point:
			 for(int i=0;i<D;i++)
			 {
				 outFile << d[line][i] << " ";
			 }
			 outFile << endl;
			 double sum = 0.0;
			 selection_sort(A, dataPoint);
			 outFile << "Knn results:\tDistance:" <<endl;
			 for (vector<int>::iterator l = A.begin();
										l != A.end();
										++l)
			 {
				 outFile << pointID[*l] << "\t\t"<< euclideandistance(&data[*l][1], dataPoint)<<endl;
				 sum += euclideandistance(&data[*l][1], dataPoint);
			 }
			 outFile<< "Total nodes searched : " <<stx::CountQueryNodes("Get")<<endl;
			 outFile << "Query Time: " << (stop - start)<< " Seconds\n";
			 outFile<< "Average k-NN distance: " << sum/k<<endl << endl;
	}

}

void Pyramid::searchquerypointsUM(double **d,int tr)
{

    pyrval pv;

    double * temp = new double[D];

    for(int i=0; i< tr; i++)
    {
    	double start = startTimer();
    	for(int j=0;j<D;j++)
    	{
    		temp[j] = d[i][j+1];
    	}
        pv=pyramidvalue_m(temp);

        outFile << "\npoint(s) found: ";
    	stx::CountQueryNodes("Start");
    	btree_type::iterator bi = bt.lower_bound(pv.pyramid + pv.height);
        while(bi != bt.end() && bi.key() <= pv.pyramid + pv.height){

        	outFile << pointID[bi.data()] << " ";

        	bi++;
        }
        outFile << endl;
        double stop = stopTimer();
        outFile<< "Total nodes searched : " <<stx::CountQueryNodes("Get")<<endl;
        outFile << "Query Time:" << (stop - start) << "  Seconds\n";
        outFile << endl;
    }

}

void Pyramid::searchqueryrangesUM(double* min, double* max, int tr)
{

	stx::CountQueryNodes("Start");
    double qr[D][2];
    double lower;
    double upper;
    for(int line=0;line<tr;line++)
	{
    	double start = startTimer();
    	stx::CountQueryNodes("Start");
		outFile << "Range Query: " << endl;
		for(int i=0; i< D; i++)
		{
			qr[i][0] = min[(line*D) + i];
			outFile << min[(line*D) + i] << ",";
			qr[i][1] = max[(line*D) + i];
			outFile << max[(line*D) + i] << "\n";

		}

	  //Use the query range to perform KNN search
		Pyramid::queryrange *q = getqueryrange_m(qr);
		outFile << "Points in range: ";
		for(int i=0; i<2*D;i++)
		{
			if(q[i].pyramid != -1)
			{
				lower = q[i].pyramid + q[i].h_low;
				upper = q[i].pyramid + q[i].h_high;
				btree_type::iterator bi = bt.lower_bound(lower);
				while(bi != bt.end() && bi.key() <= upper)
				{
					if(inside(qr, bi.data()))
					{
						outFile << pointID[bi.data()] << " ";
					}
					bi++;
				}
			}
		}
		outFile << "\n";

		outFile<< "Total nodes searched : " <<stx::CountQueryNodes("Get")<<endl;
		double stop = stopTimer();
		outFile << "Query Time: " <<  (stop - start) << " Seconds\n";
		outFile <<endl;

	}
}

void Pyramid::searchqueryKnnUM(double **d, int tr,int k)
{
	searchqueryKnnUAppM(d,tr,k);

}


void Pyramid::searchqueryranges(double* min, double* max, int tr)
{

    double qr[D][2];
    double lower;
    double upper;

    for(int line=0;line<tr;line++)
    {
    	double start = startTimer();
    	stx::CountQueryNodes("Start");
    	outFile << "Range Query: " << endl;
        for(int i=0; i< D; i++)
        {
            qr[i][0] = min[(line*D) + i];
            outFile << min[(line*D) + i] << ",";
            qr[i][1] = max[(line*D) + i];
            outFile << max[(line*D) + i] << "\n";

        }


        Pyramid::queryrange *q = getqueryrange(qr);
        outFile << "Points in range: ";
        for(int i=0; i<2*D;i++)
        {
        	if(q[i].pyramid != -1)
        	{
        		lower = q[i].pyramid + q[i].h_low;
        		upper = q[i].pyramid + q[i].h_high;
        		stx::CountQueryNodes("Start");
        		btree_type::iterator bi = bt.lower_bound(lower);
        	    while(bi != bt.end() && bi.key() <= upper)
        	    {
        	    	if(inside(qr, bi.data()))
        	    	{
        	    		outFile << pointID[bi.data()] << " ";
        	    	}
        	    	bi++;
        	    }
        	}
        }

        outFile << "\n";

        outFile<< "Total nodes searched : " <<stx::CountQueryNodes("Get")<<endl;
        double stop = stopTimer();
        outFile << "Query Time: " << (stop - start)<< " Seconds\n";
        outFile << endl;

    }
}

void Pyramid::searchqueryKnn(double **d, int tr,int k)
{


	for(int line =0; line<tr;line++)
	{
		double start = startTimer();
		 stx::CountQueryNodes("Start");
		 double dataPoint[D];
		 for(int i=0;i<D;i++)
		 {
			 dataPoint[i] = d[line][i+1]; //strip off point id
		 }

		 Pyramid::pyrval pv = pyramidvalue(dataPoint);
		 double qr[D][2];
		 double qw[D][2];
		 double tempLargestValue = 0;
		 int 	tempLargestIndex;
		 vector<int> A;
		 double dMax = 0;

		 //cout << "Before search Left " << endl;
		 searchLeft(pv, A, dataPoint, k, &dMax);
		 //cout << "Before search right " << endl;
		 searchRight(pv, A, dataPoint, k, &dMax);
		 //cout << "after search right " << endl;
		 //did we find k points in first pyrimid?
		 if( A.size() < k)
		 {
			 //make qr entire data space
			for(int i=0;i<D;i++)
			{
				qr[i][0] = 0.0;
				qr[i][1] = 1.0;
			}
		 }

		 for(int j=0; j<2*D; j++)
		 {
			 if(j != pv.pyramid)
			 {
				 //shrink qr if we have found k candidates
				 if( A.size() >= k)
				 {
					 for(int k=0;k<D;k++)
					 {

						 //calculate min in the kth dimension
						 if(dataPoint[k] - dMax < 0)
						 {
							 qr[k][0] = 0;
						 }
						 else
						 {
							 qr[k][0] = dataPoint[k] - dMax;
						 }

						 //calculate max in the kth dimension
						 if(dataPoint[k] + dMax > 1)
						 {
							 qr[k][1] = 1;
						 }
						 else
						 {
							 qr[k][1] = dataPoint[k] + dMax;
						 }
					 }
				 }


				 //find intersection of qr with j
				 double hlow, hhigh;
				 for(int i=0; i<D; i++)
				 {
					 qw[i][0] = qr[i][0] - 0.5;
					 qw[i][1] = qr[i][1] - 0.5;
				 }
				 if(intersect(qw,j,&hlow,&hhigh))
				 {
					 //update A and dMax
					double lower = j + hlow;
					double upper = j + hhigh;
					btree_type::iterator bi = bt.lower_bound(lower);
					while(bi != bt.end() && bi.key() <= upper)
					{
						//check if bi.data() needs to be added to A
						if(A.size() < k)
						{
							 //add point to A then check if we need to update dMax
							 //cout << "adding " << bi.data() << " because for A.size: " << A.size();
							 A.push_back(bi.data());

							 if(euclideandistance(&data[bi.data()][1], dataPoint) > dMax )
							 {
								 dMax = euclideandistance(&data[bi.data()][1], dataPoint);
							 }
						}
						else
						{
							//if distance between bi  and q < dMax update A and dMax
							 if(euclideandistance(&data[bi.data()][1], dataPoint) < dMax && !in(A, bi.data()))
							 {
								 //update A
								 tempLargestValue = 0;
								 for(int i=0;i<A.size(); i++)
								 {
									 if(euclideandistance(&data[A[i]][1], dataPoint) > tempLargestValue)
									 {
										 tempLargestValue = euclideandistance(&data[A[i]][1], dataPoint);
										 tempLargestIndex = i;
									 }
								 }

								 A[tempLargestIndex] = bi.data();

								 //update dMax
								 tempLargestValue = 0;
								 for(int i=0;i<A.size(); i++)
								 {
									 if(euclideandistance(&data[A[i]][1], dataPoint) > tempLargestValue)
									 {
										 tempLargestValue = euclideandistance(&data[A[i]][1], dataPoint);
									 }
								 }
								 dMax = tempLargestValue;
							 }
						}

						bi++;
					}
				 }
			 }
		 }
		 outFile << "k-NN Query " << endl;
		 //print out query point:
		 for(int i=0;i<D;i++)
		 {
			 outFile << dataPoint[i] << " ";
		 }
		 outFile << endl;
		 double sum = 0.0;
		 outFile << "Knn results:\tDistance:" <<endl;
		 double stop = stopTimer();
		 selection_sort(A, dataPoint);
		 for (vector<int>::iterator l = A.begin();
									l != A.end();
									++l)
		 {
			 outFile << pointID[*l] << "\t\t"<< euclideandistance(&data[*l][1], dataPoint)<<endl;
			 sum += euclideandistance(&data[*l][1], dataPoint);
		 }
		 outFile<< "Total nodes searched : " <<stx::CountQueryNodes("Get")<<endl;
		 outFile << "Query Time: " << (stop - start) <<  " Seconds\n";
		 outFile<< "Average k-NN distance: " << sum/k<<endl;
		 outFile << endl;
	}

}

void Pyramid::searchLeft(pyrval pv, vector<int>& A, double q[], int k , double* dMax)
{

	 btree_type::iterator bi = bt.upper_bound(pv.pyramid + pv.height);
	 bi--;
	 double tempLargestValue = 0;
	 int 	tempLargestIndex;
	 while( bi.key() >= pv.pyramid )
	 {
		  if(A.size() < k && bi.key() >= pv.pyramid && !in(A, bi.data()))
		 {
			 //add point to A then check if we need to update dMax
			 A.push_back(bi.data());

			 if(euclideandistance(&data[bi.data()][1], q) > *dMax )
			 {
				 *dMax = euclideandistance(&data[bi.data()][1], q);
			 }
		 }
		 else if(bi.key() >= pv.pyramid)
		 {
			 //if distance between bi  and q < dMax update A and dMax
			 if(euclideandistance(&data[bi.data()][1], q) < *dMax && !in(A, bi.data()))
			 {
				 //update A
				 tempLargestValue = 0;
				 for(int i=0;i<A.size(); i++)
				 {
					 if(euclideandistance(&data[A[i]][1], q) > tempLargestValue)
					 {
						 tempLargestValue = euclideandistance(&data[A[i]][1], q);
						 tempLargestIndex = i;
					 }
				 }
				 A[tempLargestIndex] = bi.data();

				 //update dMax
				 tempLargestValue = 0;
				 for(int i=0;i<A.size(); i++)
				 {
					 if(euclideandistance(&data[A[i]][1], q) > tempLargestValue)
					 {
						 tempLargestValue = euclideandistance(&data[A[i]][1], q);
					 }
				 }
				 *dMax = tempLargestValue;

			 }

		 }
		 if(bi == bt.begin())
		 {
			 break;
		 }
		 bi--;
	 }

}

void Pyramid::searchRight(pyrval pv, vector<int>& A, double q[], int k, double* dMax)
{

	 btree_type::iterator bi = bt.lower_bound(pv.pyramid + pv.height);
	 double tempLargestValue = 0;
	 int 	tempLargestIndex;
	 while( bi.key() < pv.pyramid +1.0)
	 {
		 if(A.size() < k && bi.key() <= pv.pyramid +1 && !in(A, bi.data()))
		 {
			 //add point to A then check if we need to update dMax
			 A.push_back(bi.data());

			 if(euclideandistance(&data[bi.data()][1], q) > *dMax )
			 {

				 *dMax = euclideandistance(&data[bi.data()][1], q);
			 }
		 }
		 else if(bi.key() < pv.pyramid +1.0)
		 {

			 //if distance between bi  and q < dMax update A and dMax
			 if(euclideandistance(&data[bi.data()][1], q) < *dMax && !in(A, bi.data()))
			 {
				 //update A
				 tempLargestValue = 0;
				 for(int i=0;i<A.size(); i++)
				 {
					 if(euclideandistance(&data[A[i]][1], q) > tempLargestValue)
					 {
						 tempLargestValue = euclideandistance(&data[A[i]][1], q);
						 tempLargestIndex = i;
					 }
				 }
				 A[tempLargestIndex] = bi.data();

				 //update dMax
				 tempLargestValue = 0;
				 for(int i=0;i<A.size(); i++)
				 {
					 if(euclideandistance(&data[A[i]][1], q) > tempLargestValue)
					 {
						 tempLargestValue = euclideandistance(&data[A[i]][1], q);
					 }
				 }
				 *dMax = tempLargestValue;
			 }

		 }

		 bi++;
		 if(bi == bt.end())
		 {
			 break;
		 }
	 }
}
//checks to see if y is in vector x
bool Pyramid::in(vector<int> x, int y)
{
	for(vector<int>::iterator i = x.begin(); i != x.end(); ++i)
	{
		if(*i == y)
			return true;
	}

	return false;
}

void Pyramid::searchquerypointsSS(double **d, int tr)
{
	double * temp = new double[D];

    for(int i=0; i< tr; i++)
    {
    	double start = startTimer();
    	for(int j=0;j<D;j++)
    	{
    		temp[j] = d[i][j+1]; // strip off point id
    	}

    	//find temp in data
    	bool found = true;
    	outFile << "Points found: ";
    	for(int j=0;j<totalrecords;j++)
    	{
    		found = true;
    		for(int k=0;k<D;k++)
    		{
    			if( data[j][k+1] != temp[k])
    			{
    				found = false;
    				break;
    			}
    		}
    		if(found)
    		{
    			outFile << pointID[j] << " ";
    		}
    	}
    	outFile << endl;
    	outFile<< "Total nodes searched : " <<totalrecords<<endl;
    	double stop = stopTimer();
    	outFile << "Query Time: " << (stop - start) <<" Seconds\n";
        outFile << endl;
    }
}

void Pyramid::searchqueryrangesSS(double *min, double *max, int tr)
{
    double qr[D][2];


    for(int line=0;line<tr;line++)
    {
    	double start = startTimer();
        outFile << "Range Query: " << endl;
        for(int i=0; i< D; i++)
        {
            qr[i][0] = min[(line*D) + i];
            outFile << min[(line*D) + i] << ",";
            qr[i][1] = max[(line*D) + i];
            outFile << max[(line*D) + i] << "\n";
        }


    	bool found = true;
    	outFile << "Points found: ";
    	for(int j=0;j<totalrecords;j++)
    	{
    		found = true;
    		for(int k=0;k<D;k++)
    		{
    			if( data[j][k+1] < qr[k][0] || data[j][k+1] > qr[k][1])
    			{
    				found = false;
    				break;
    			}
    		}
    		if(found)
    		{
    			outFile << pointID[j] << " ";
    		}
    	}

    	outFile << "\n";

    	outFile<< "Total nodes searched : " <<totalrecords<<endl;
        double stop = stopTimer();
        outFile << "Query Time: "  <<  (stop - start) << " Seconds\n";
        outFile << endl;
    }
}

void Pyramid::searchqueryKnnSS(double **d, int tr, int k)
{
	for(int line =0; line<tr;line++)
	{

		 double sum = 0;
		 double start = startTimer();

		 double dataPoint[D];
		 for(int i=0;i<D;i++)
		 {
			 dataPoint[i] = d[line][i+1]; //strip off point id
		 }

		 double tempLargestValue = 0;
		 int 	tempLargestIndex;
		 vector<int> A;
		 double dMax = 0;

		 for(int i=0;i<totalrecords;i++)
		 {

				 if(A.size() < k)
				 {//add this point;

					 A.push_back(i);

					 //check to see if we should update dMax;
					 if(euclideandistance(dataPoint, &data[i][1]) > dMax)
					 {
						 dMax = euclideandistance(dataPoint, &data[i][1]);
					 }
				 }
				 else
				{
					//if distance between bi  and q < dMax update A and dMax
					 if(euclideandistance(dataPoint, &data[i][1]) < dMax )
					 {
						 //update A
						 tempLargestValue = 0;
						 for(int k=0;k<A.size(); k++)
						 {
							 if(euclideandistance(&data[A[k]][1], dataPoint) > tempLargestValue)
							 {
								 tempLargestValue = euclideandistance(&data[A[k]][1], dataPoint);
								 tempLargestIndex = k;
							 }
						 }

						A[tempLargestIndex] = i;

						 //update dMax
						 tempLargestValue = 0;
						 for(int k=0;k<A.size(); k++)
						 {
							 if(euclideandistance(&data[A[k]][1], dataPoint) > tempLargestValue)
							 {
								 tempLargestValue = euclideandistance(&data[A[k]][1], dataPoint);
							 }
						 }
						 dMax = tempLargestValue;
					 }
				}
		 }


		 outFile << "Knn results:\tDistance:" <<endl;
		 double stop = stopTimer();
		 //sort A
		 selection_sort(A, dataPoint);
		 for (vector<int>::iterator l = A.begin();
									l != A.end();
									++l)
		 {
			 outFile << pointID[*l] << "\t\t"<< euclideandistance(&data[*l][1], dataPoint)<<endl;
			 sum += euclideandistance(&data[*l][1], dataPoint);
		 }

		 outFile<< "Total nodes searched : " <<totalrecords<<endl;

		 outFile << "Query Time: " << (stop - start)<< " Seconds\n";
		 outFile<< "Average k-NN distance: " << sum/k<<endl;
		 outFile << endl;
	}
}

Pyramid::pyrval Pyramid::pyramidvalue(double values[]){
    int dmax=0;
	    int pyramidnumber=0;
	    double height=fabs(0.5-values[0]);
	    for(int j=1; j<D; j++){
	        if(height<fabs(0.5-values[j]))
	        {
	            dmax = j;
	            height = fabs(0.5-values[j]);
	        }
	    }
	    if(values[dmax] < 0.5)
	    {
	        pyramidnumber = dmax;
	    }
	    else
	    {
	        pyramidnumber = D + dmax;
	    }

	    Pyramid::pyrval pv(pyramidnumber,height);
	    return pv;
    
}

//A function to calculated pyramid value
Pyramid::pyrval Pyramid::pyramidvalue(std::vector<double> values){

	    int dmax=0;
	    int pyramidnumber=0;
	    double height=fabs(0.5-values[0]);
	    for(int j=1; j<D; j++){
	        if(height<fabs(0.5-values[j]))
	        {
	            dmax = j;
	            height = fabs(0.5-values[j]);
	        }
	    }
	    if(values[dmax] < 0.5)
	    {
	        pyramidnumber = dmax;
	    }
	    else
	    {
	        pyramidnumber = D + dmax;
	    }

	    Pyramid::pyrval pv(pyramidnumber,height);
	    return pv;

	}


//A function to calculate the pyramid number (the actual pyramid number and not the main dimension of the pyramid) and height low and height high for range search
Pyramid::queryrange* Pyramid::getqueryrange(double qr[][2]){

    Pyramid::queryrange *q= (queryrange *)malloc(sizeof(queryrange)*2*D);

    double qw[D][2];

    //pyramid query
    for(int i=0; i<D; i++)
    {
    	qw[i][0] = qr[i][0] - 0.5;
    	qw[i][1] = qr[i][1] - 0.5;
    }


	for(int i=0;i<2*D;i++)
	{

		double hlow, hhigh;

		if(intersect(qw,i,&hlow,&hhigh))
		{
						q[i].pyramid = i;
						q[i].h_high = hhigh;
						q[i].h_low = hlow;

		}
		else
		{
			q[i].pyramid = -1;
		}
	}
    return q;
}

//A function to calculate absolute min of given values
double Pyramid::MIN(double *r)
{
	if(r[0]*r[1]<=0) return 0;
	else return fabs(r[0])<fabs(r[1])? fabs(r[0]):fabs(r[1]);
}

//A function to calculate absolute max of given values
double Pyramid::MAX(double *r)
{
	return fabs(r[0])>fabs(r[1])? fabs(r[0]):fabs(r[1]);
}



//A function to check whether pyramid i, intersects with qw[][]. If it intersects it return true with values of height low and height high in hlow and hhigh.
//This function uses DetermingRange internally to calculate hlow and hhigh
//Adopted from P+ tree source code
bool Pyramid::intersect(double qw[][2], int i, double* hlow, double* hhigh){


    int j;
	if(i<D)
	{
		for(j=0;j<D && qw[i][0]<=-MIN(qw[j]);j++);
		if(j==D)
		{
			double qb[D][2];
			for(int k=0;k<D;k++)
			{
				for(int l=0;l<2;l++)
				{
					qb[k][l] = qw[k][l];
				}
			}
			if (qb[i][1]>0) qb[i][1]=0;
			DeterminRange(qb,i,hlow,hhigh);
			return 1;
		}
		else return 0;
	}

	else
	{

		for(j=0;j<D && qw[i-D][1]>=MIN(qw[j]);j++);
		if(j==D)
		{
			double qb[D][2];
			for(int k=0;k<D;k++)
			{
				for(int l=0;l<2;l++)
				{
					qb[k][l] = qw[k][l];
				}
			}
			if (qb[i-D][0]<0) qb[i-D][0]=0;

			DeterminRange(qb,i-D,hlow,hhigh);
			return 1;
		}
		else return 0;
	}
    return 0;
}

//A function to height low (hlow) and height high (hhigh) in the intersection ofpyramid i and qw[][]
//Adopted from P+ tree source code
void Pyramid::DeterminRange(double qw[][2], int i, double* hlow, double* hhigh){


	int j,k;
	double qtemp;
	*hhigh=MAX(qw[i]);
	for(j=0;j<D && qw[j][0]<=0 && 0<=qw[j][1]; j++);
	if(j==D) *hlow=0;
	else
	{
		*hlow=0.0;
		for(k=0;k<D;k++)
		{
			qtemp=MIN(qw[k])>MIN(qw[i])? MIN(qw[k]):MIN(qw[i]);
			if(qtemp > *hlow) *hlow=qtemp;
		}
	}
}

//A function to find approximate medians using histograms. This method calls functions of ApproximateMedian class
void Pyramid::hbam(int tr, int d, double **val){
    am.initialize(tr,d,val);
    am.findapproximatemedians();
    
}

//A function to find approximate medians using algorithm. This method calls functions of ApproximateMedian class
void Pyramid::algam(int tr, int d, double **val){
    am.initialize(tr,d,val);
    am.findapproximatemediansusingAMN();
}

void Pyramid::findmedian(int tr, int d, double **val){
    am.initialize(tr,d,val);
    am.findmedian();
}


//A function to calculated pyramid value
//using approximate median
Pyramid::pyrval Pyramid::pyramidvalue_am(double values[]){

        double * convertedval;
        
        convertedval= am.convertdatapoint(values);


        
	    int dmax=0;
	    int pyramidnumber=0;
	    double height=fabs(0.5-convertedval[0]);
	    for(int j=1; j<D; j++){
	        if(height<fabs(0.5-convertedval[j]))
	        {
	            dmax = j;
	            height = fabs(0.5-convertedval[j]);
	        }
	    }
	    if(convertedval[dmax] < 0.5)
	    {
	        pyramidnumber = dmax;
	    }
	    else
	    {
	        pyramidnumber = D + dmax;
	    }

	    Pyramid::pyrval pv(pyramidnumber,height);
	    return pv;

	}
//A function to calculated pyramid value
//using approximate median
Pyramid::pyrval Pyramid::pyramidvalue_m(double values[]){

		double * convertedval;

		convertedval= am.convertdatapoint(values);

		int dmax=0;
		int pyramidnumber=0;
		double height=fabs(0.5-convertedval[0]);
		for(int j=1; j<D; j++){
			if(height<fabs(0.5-convertedval[j]))
			{
				dmax = j;
				height = fabs(0.5-convertedval[j]);
			}
		}
		if(convertedval[dmax] < 0.5)
		{
			pyramidnumber = dmax;
		}
		else
		{
			pyramidnumber = D + dmax;
		}

		Pyramid::pyrval pv(pyramidnumber,height);
		return pv;

	}


//A function to calculate the pyramid number (the actual pyramid number and not the main dimension of the pyramid) and height low and height high for range search,
//using the approximate median
Pyramid::queryrange* Pyramid::getqueryrange_am(double qval[][2]){


    Pyramid::queryrange *q= new queryrange[2*D];

    double qw[D][2];
    double** qr=    am.convertqueryrange(qval);

    //pyramid query
    for(int i=0; i<D; i++)
    {
    	qw[i][0] = qr[i][0] - 0.5;
    	qw[i][1] = qr[i][1] - 0.5;
    }

	for(int i=0;i<2*D;i++)
	{

		double hlow, hhigh;

		if(intersect(qw,i,&hlow,&hhigh))
		{
				q[i].pyramid = i;
				q[i].h_high = hhigh;
				q[i].h_low = hlow;

		}
		else
		{
			q[i].pyramid = -1;
		}
	}

    return q;
}

Pyramid::queryrange* Pyramid::getqueryrange_m(double qval[][2]){


    Pyramid::queryrange *q= new queryrange[2*D];

    double qw[D][2];

    double** qr = am.convertqueryrange(qval);

    //pyramid query
    for(int i=0; i<D; i++)
    {
    	qw[i][0] = qr[i][0] - 0.5;
    	qw[i][1] = qr[i][1] - 0.5;
    }


	for(int i=0;i<2*D;i++)
	{

		double hlow, hhigh;

		if(intersect(qw,i,&hlow,&hhigh))
		{

			q[i].pyramid = i;
			q[i].h_high = hhigh;
			q[i].h_low = hlow;
		}
		else
		{
			q[i].pyramid = -1;
		}
	}

    return q;
}

//checks that data point id is inside the query range qr
bool Pyramid::inside(double qr[][2], int id){
	for(int i=0;i<D;i++)
	{
		if(qr[i][0] > data[id][i+1])
		{
			return false;
		}
		if(qr[i][1] < data[id][i+1])
		{
			return false;
		}
	}
	return true;
}

//Calculate and return euclidean distance between two points
double Pyramid::euclideandistance(double a[], double b[]){


        double sum=0.0;
        for(int i=0; i<D; i++)
        {
            sum = sum + pow( (a[i] - b[i]) , 2 );
        }
        return  sqrt(sum);
    }

//This function sorts Size elements of the array A located at positions Left, Left+Step,
//Left+2*Step.....
void Pyramid::selection_sort(vector<int>& A, double point[]){
    int min=0;
    int swap1 = 0, swap2 = 0;
    for(int i=0; i<A.size(); i++)
    {
        min = i;
        for(int j= i +1; j< A.size(); j++)
        {
        	euclideandistance(&data[ A[j] ][1], point);
            if(euclideandistance(&data[ A[j] ][1], point) < euclideandistance(&data[A[min]][1], point)){
                min = j;
            }
        }
        swap1 = A[i];
        swap2 = A[min];
        A[i]=swap2;
        A[min]=swap1;
    }
}


